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L POLYMER DYNAMICS AND FLOW PROPERTIES WE WANT TO UNDERSTAND: MOTIVATION AND GOALS 

o : 

f*^ A. Challenges in polymer dynamics under flow 

(N '■ 

J> , Polymer molecules differ from simple fluids in several aspects: they are extremely diverse in structure (they can have a linear, 
branched, ring-like, or block copolymer structure), they can be characterized by a molecular weight distribution, and they are 
capable of exhibiting a huge number of configurations implying that a large number of degrees of freedom should be accounted 
for in any molecular modeling approach. As a result, polymers exhibit properties which are totally distinct from those of the 
simpler Newtonian liquids. The drag reduction phenomenon (the substantial reduction in pressure drop during the turbulent 
,— — | . flow of a Newtonian liquid when a very small amount of a flexible polymer is added), their unique rheological properties 
(shear thinning and normal stress differences in simple shear, strain hardening in elongation, complex viscosity, anisotropy in 
thermal conductivity and diffusivity), and a plethora of other phenomena associated with their elastic character are only a few 
manifestations of the departure of their behavior from the Newtonian one QUI]. Of particular importance from a mechanical 
or fluid dynamics point of view is their viscoelasticity quantifying the irreversible conversion of the work needed for their 
deformation to heat loss but also their capability to store part of this work as elastic energy. It is a property closely related to 
the multiplicity of time and length scales characterizing dynamics and structure in these fluids. Thus, even in the viscous regime 
(Wi -C 1, where Wi is the Weissenberg number empirically defined as Wi = r p j with r p being the longest relaxation time and 
7 the flow rate), the flow can still be strong enough for several degrees of freedom not to be close to equilibrium giving rise to 
interesting rheological properties also there [2], especially for high molecular weight polymers. 

Understanding relaxation processes and structure development occurring over these multiple scales is a prerequisite for deriv- 
ing reliable constitutive equations connecting the stresses developing in these materials in terms of the imposed flow kinematics 
and certain molecular parameters or functions, and for computing polymer flows 0, H, 0, 0, S, HI [Toll . It is only through 
a comprehensive understanding across scales that one can hope to build the relationship between polymer molecular structure, 
conformation and architecture and macroscopic rheological response. In addition to experiments and theories, molecular simula- 
tions can play a significant role by providing high resolution calculations especially on the crossover from small-to-intermediate 
scales. This chapter is devoted to a brief discussion of some of the emerging multi-scale simulation approaches in the recent 
research literature on nonequilibrium systems, with emphasis on those based on well-founded theoretical frameworks. Our goal 
is to demonstrate that, with the help of and guidance from recent advances in the field of nonequilibrium statistical mechanics 
and thermodynamics, this highly demanding endeavor (modeling across scales) can lead to simulation methodologies that have 
been elevated from simple, brute-force computational experiments to systematic tools for extracting complete, redundancy-free, 
and consistent coarse-grained information for the flow dynamics of polymeric systems 1 1 111 . 
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5_] ■ B. Modeling polymer dynamics beyond equilibrium 

a _ 

Describing macromolecular configurations under nonequilibrium conditions is an extremely difficult problem, which usually 
requires simplified models for analytical or numerical studies [2]; such simple models have contributed enormously to our 
understanding of polymer rheology and mechanics. For a review on proposed models, see e .g. fl^fl^fl^fl^flql . while for a 
review on available simulation tools addressing different time and length scales, see ifl^fTllHIHMrH |23ll. Fig. Q] shows 
schematic of the pertinent models for polymer solutions and melts depending on the length- and time-scale of interest. 

In general, examples of macroscopic constitutive equations employed to calculate polymer flow behavior include e.g. confor- 
mation tensor models like the Giesekus, Maxwell and FENE models but also the more recently proposed Pom-Pom, CCR and 
Rolie-Poly models lfl5l[l6ll24ll25ll26ll27ll . For an overview, see (HIU as well as the contributions by A.N. Beris and E. Mitsoulis 
in this volume. Most of these constitutive equations have been derived (or inspired) by simple mechanical models of polymer 
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FIG. 1: Different models of polymer dynamics are shown schematically for solutions and melts. 



motion. They are mesoscale, kinetic theory models (based, e.g., on the dumbbell, FENE, bead-spring chain, and bead-rod chain 
analogues) capable of accounting for some important aspects of polymer dynamics, like chain stretching, non-affine deforma- 
tion, diffusion, and hydrodynamic drag forces, either separately or all together ll28tl . In polymer solutions, in particular, one 
should account for the solvent-mediated effect between beads (known as hydrodynamic interaction), which is usually modeled 
with the hydrodynamic resistance matrix |2j,|29J]. Recently, an efficient simulation of solvent dynamics with the help of either 
the lattice Boltzmann lf3Qt l3lll or the stochastic rotation dynamics J32I l33ll methods has been proposed with a suitable coupling 
of the bead dynamics. Besides their efficiency, these combined methods also incorporate non-uniformities and fluctuations in 
internal flow fields. 

Usually, the mesoscopic, kinetic models are considered to be well-suited for predicting dynamical properties of polymer 
solutions on macroscopic scales. Details of the fast solvent dynamics are in most cases irrelevant for macroscopic properties. 
Exceptions are polyelectrolytes, where the motion of counterions in the solvent can have a major influence on polymer conforma- 
tion. Therefore, more microscopic models of polyelectrolytes with explicit counterions are sometimes employed IB4I1 (see also 
the contribution by M. Muthukumar in this volume). Another exception is the dynamics of individual biopolymers, e.g. protein 
folding, which is modeled with an all atomistic model including an explicit treatment of the (water) solvent molecules 113511 - 

For typical polymer melts, dumbbell models are inappropriate since they fail to account for the essential role of entanglements 
on the long-time dynamics. Successful mesoscopic, mean-field descriptions of entangled melts are offered by the reptation 
models 11411 . Several modifications of the original Doi-Edwards-de Gennes reptation theory have been prop osed in the last 
years in order to improve comparison to experiments in the linear and nonlinear flow regimes |36[ l37i f38l l39l l40l l4ll r42l l43ll . 
Recently, slip-link models have also been proposed Ir44lr45lr4q, l47. 48, 49, 50], providing a slightly more detailed description of 
entanglements, and which agree well with available experimental results. 

At the microscopic level, polymer melts are modeled as multi-chain systems, see e.g. J2, EH EH Hill • For example, all-atom 
or united atom force fields, accounting explicitly for bond angle bending and torsion angle contributions (in addition to bond 
stretching and intermolecular interactions) 1521 15311 . are available. Different united-atom force fields are reviewed and compared 
e.g. in I5lll54ll55ll . From such detailed atomistic molecular dynamics (MD) simulations, the linear viscoelastic properties can be 
computed by Green-Kubo relationships lf56l l57l l58ll . Also less detailed bead-spring models are available; a prototypical model 
is that of Kremer and Grest ll59ll (or variants thereof l60|p . which neglects chemical details and instead focuses on the interplay 
between chain connectivity and excluded volume effects. The non-linear regime can be studied by nonequilibrium molecular 
dynamics (NEMD) simulations, that directly address flow effect on polymer structure and conformation, both in shear and 
planar elongation (see e.g. lfl2ll6^loul^lo^l^lo^loalo7n and references therein). They are based on flow-adapted boundary 
conditions such as those proposed by Lees-Edwards for planar shear [68] and by Kraynik-Reinelt for planar elongation j&J. 



C. Challenges in standard simulations of polymers in flow 

Despite the enormous advances in the field of molecular simulations ffl E^l2lll2^l^l7lll^l73Tl . predicting the macroscopic 
flow properties of polymers from their underlying microstructure presents still major challenges 17411 . Available MD and NEMD 
algorithms can address only time scales on the order of a few microseconds at most, implying that only short up to moderately 
entangled polymers can be studied in full atomistic detail in a brute-force manner. Extending the simulations to longer, truly 
entangled polymers is a first big task. Extending these flows to mixed or inhomogeneous flows is another big challenge. Among 
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others, such a development would help understand the origin of interfacial slip and its mode (localized slip versus global slip) 
in the flow past a solid substrate. On the other hand, with the introduction of the revolutionary set of chain connectivity altering 
moves, extremely powerful Monte Carlo algorithms have been developed which have helped overcome the issue of the thermal 
equilibration of long polymers even at beyond equilibrium conditions IT751 l76i PhI p78i p79tl . With the help of the end-bridging 
and double bridging moves, for example, truly long polyethylene (linear and branched) and polybutadiene systems have been 
equilibrated over the last years, which opened also the way to their topological analysis for the identification of entanglements 

m mm mm m m. 

Arguably the biggest challenge in polymer simulation under nonequilibrium conditions is to build well-founded multi-scale 
tools which can bridge the gap between microscopic information and macroscopically manifested viscoelastic properties, prefer- 
ably through a constitutive equation founded on the microscopic mode l fill 17011 . Simulations of metals face the same problem 
where again the objective is concurrent length scale simulations l87i l88ll ; to some extent, it is also relevant to simple fluids 
118911 - For polymers, additional motivation stems from the increased interest in polymer mixtures and interfaces 19011 resulting in 
morphology development at the nanoscale. 

Here, we aim to shortly review some recent coarse-graining and multi-scale methods (see I23I HH |79l |9jl H3, HH l95ll 
and also [96, ^Ulsil] for recent reviews of such methods for polymers and in a more general context, respectively), but also 
to put forward some new ideas for addressing such issues, which could eventually allow to model the macroscale quantities 
of interest by a suitable coarse-graining procedure. We will see that, if one is guided by nonequilibrium statistical mechanics 
and thermodynamics, it is possible to design well-founded multi-scale modeling tools that can link microscopic models with 
macroscopic constitutive equations. Such multi-scale modeling t ools benef it from recently proposed approaches for static coarse 
graining which are mainly built on potentials of mean force lT79L 1 1 OOL 1 1 ill - For dynamical properties, however, coarse-grained 
models need to account for dissipative effects that arise due to fast d egrees of free dom that are eliminated (for example via 
projection operators) in favor of the remaining slowly varying ones ll56ll58l[T02ill03ll . In a flow situation, these slow dynamical 
variables depart from their values in the quiescent fluid, while all other (faster) degrees of freedom track the evolution of 
the structural parameters; i.e., they are assumed to be in local equilibrium subject to the constraints imposed by the values 
of the structural parameters at all times. The proper definition of the set of the state variables, effectively representing the 
nonequilibrium states is the first key to the success of such an approach. Linking the microscopic model with a macroscopic 
model built on these slowly relaxing variables is the second key; as we will discuss in the next sections of this chapter, this is best 
addressed by getting guidance from a nonequilibrium statistical thermodynamics framework proposing a fundamental evolution 
equation for the macroscopic model in terms of the chosen structural (dynamical) variables. 

II. COARSE-GRAINED VARIABLES AND MODELS 

We start with a microscopic polymer model, whose state is specified by a point in 6./V-dimensional phase space, z 6 T with 
z = (ri, . . . , rjy; pi, . . . , Pat) a short notation for the positions and momenta of all N particles. The model is described by the 
microscopic Hamiltonian H (z) with inter- and intra-molecular interactions. The coarse-grained model eliminates some of the 
(huge number of) microscopic degrees of freedom. The level of detail that is retained is specified by the choice of coarse-grained 
variables x = (sci, . . . , x ncG ) with 



Xk = (life) = J dzll fe (z)p(z), k=l, ...,ucg, (1) 

where p(z) denotes the probability distribution on V and the phase space functions IT^z) are the instantaneous values of the 
coarse-grained variables in the microstate z. 

Instead of the full, microscopic distribution p(z), the coarse-grained model is already specified by the reduced probability 
distribution p(x) = (<5(x — n(z))). Knowledge of p(x) allows to calculate averages of quantities a(n(z)) via (a(n)) = 
/ r dz a(n(z))p(z) = J dxa(x)p(x). Instead of p(x), coarse-grained models are often described by the so-called potentials 
of mean force U m i(x) = — fceTlnp(x) which formally replace the Hamiltonian in the calculation of equilibrium, canonical 
averages. However, U m { is an effective free energy difference and therefore depends on the thermodynamic state of the system. 
It contains in general effective many-body in teractions that arise by partially integrating out microscopic degrees of freedom. 
For recent reviews see e.g. JHCMIIolEol. 



Different sets of coarse-grained variables have been suggested in the literature and are briefly reviewed in the following. 
Usually, one is interested in a drastic reduction of microscopic complexity, so hqg *C 6N. However sometimes, e.g. for 
equilibrating atomistic systems, a coarse-grained model with a modest reduction (ncG only a factor 5 to 10 smaller than 6 AO 
might be useful. We emphasize that the mapping (fTJ as well as the corresponding probability p(x) is not restricted to equilibrium 
situations. We proceed with a short review of coarse-grained variables and models that capture different levels of detail and 
briefly discuss their static and dynamic properties. 
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A. Beads and superatoms 

Coarse graining to the level of superatoms CM EM EM EM EH EH EH EH Ell EH EM EH EH EH EU 
is a method followed when one wishes to reduce chemical complexity in a polymer chain without loosing the chemical iden- 
tity of the molecule. According to this, a certain number of atoms or repeat units along the chain are grouped together into 
"superatoms" connected by effective bonds and governed by softer or smoother effective non-bonded interactions. The result- 
ing (usually linear) chain sequence is simpler and amenable to fast thermal equilibration through application of state-of-the-art 
Monte Carlo algorithms slightly modified to account for the presence of the few different species along the chain. The method 
has been widely applied to reduce complexity and permit the molecular simulation of a number of polymers. Typical examples 
include polystyrene, poly(ethylene terepthalate), polycarbonates, polyphenylene dendrimers, even DNA. It involves computing 
the effective intra and intermolecular potentials among superatoms such that the coarse-grained model reproduces as faithfully as 
possible the structural, configurational and thermodynamic properties of the original atomistic polymer model. For vinyl chains 
presenting sequences of methylene (CH2) and pseudo-asymmetric methyne (-CHR) groups along their backbone, the method 
should also account for the isotactic, syndiotactric or atactic stereochemistry of the polymer, based on the succession of meso 
(m) and racemo (r) diads (see Fig. [2]). 

According to Zwicker and Lovett lll2lll . if all interactions with potentials V^ n ' (ri, r 2 , . . . , r/v), with denoting the position 
vector of the i-th particle, in an iV-atom molecular system consist of rt-body and lower terms, then the system can be completely 
described by the knowledge of all 71-order correlation functions jW (n, r%, ■ . ■ , r^r) and lower. Since, in practice, the complete 
determination of the n-point correlation functions is a huge task for N > 2 and n > 4, the calculations are usually limited to 
correlation functions that depend only on a single coordinate. For polymers where potential functions are usually separated 
into intra- and intermolecular ones, examples of such correlation functions include typically the radial distribution function, 
the distribution of bond lengths, the distribution of bending angles, and the distribution of dihedral angles. The coarse-grained 
potential then should be chosen such that it matches the distributions of all possible bond lengths, bond angles, and torsional 
angles, and of all intermolecular pair distribution functions, as extracted from simulations with the corresponding atomistic 
model. 

For a distribution function that depends on a single co ordinate, t he co rresponding effective potential can be computed through 
a method called iterative Boltzmann inversion method 

CM EM EM 1 1 2211 . aimed at matching the distribution of the relevant 
degrees of freedom (called target distribution) between the chosen coarse-grained model and the initial atomistic model for the 
polymer; the latter are usually extracted from accurate, brute-force MD simulations on short homologues. The method uses 
the differences in the potentials of mean force between the distribution function generated from a guessed potential and the 
true distribution function to improve the effective potential ite rativ ely. Qualitative arguments for the conditions under which 
convergence should be expected have been discussed by Soper lll22ll . 

Naive use of the coarse- grained potent ials in standard molecular dynamics simulations leads to wrong predictions of diffusion 
and relaxation processes lU04l Il23l Il24ll . A simple, empirical method for relating the dynamics of super atom s at the coarse- 
grained level with the dynamics of true atomistic units at the atomic level uses a time rescaling factor Il25l Il26ll . Within 
this method, the effective potentials are used in the reversible equations of motion of classical mechanics to perform standard 
molecular dynamics simulations and then an effort is made to match the mean-square displa cements o f the relevant structural 
units in the atomistic and coarse-grained models, both in amplitude and slope. Noid et al. lll27l Il28ll formulate consistency 
criteria that should be obeyed when using coarse-grained potentials in equilibrium dynamical simulations. 

Coarse graining to the level of superatoms has drawn a lot of attention in the last years mainly because of the capability to 
account for the correct stereochemical sequence of the repeat units. Despite the success of effective pair potentials in reproducing 
many of the structural properties of the corresponding atomistic system, however, their use in actual simulations is accompanied 
by a number of thermodynamic inconsistencies: (a) They perform well for the particular physical properties for which they were 
developed. For example, the value of pressure as computed by using the virial theorem from the effective potentials optimized 
with the iterative Boltzmann method is higher than what is observed in the atomistic system, unless an attractive perturbation 
potential is added (ramp correction to the pressure) and the potential is re-optimized. Given that in integral equation methods 
the mechanical properties of a syste m (such as pre ssure, energy and compressibility) are fixed by the singlet and pair number 
densities along with proper closures lll2l[|l29lll30ll . such a pressur e inco nsistency should be related to the degree of sensitiv ity of 
the site-site pair correlation function to the effective pair potential Jl22ll . This is in line with the simulations of Jain et al. [131] 
who showed that, although there is a one-to-one correspondence between the structure of a liquid (i.e., the pair correlation 
function) and its pair- wise additive intermolecular potential (Henderson's theorem II 12911 ). the convergence of potentials obtained 
by standard inversion procedures is extremely slow: although the repulsive part of the potential converges rapidly, its attractive 
part (to which, e.g., the internal energy and pressure are primarily sensitive) converges slowly, (b) Effective potentials are in 
general not transferable; they are state-point (e.g., temperature and pressure) dependent. In same cas es, it has been noticed that 
temperature changes at about the same density do not drastically affect their paramete rizati on lll32lll33ll . A newly developed 
effective force coarse-graining seems to improve transferability to other state points lll34ll . Developing fully self-consistent 
and transferable potentials at any arbitrary level of coarse graining remains still a challenge, (c) The coarse-grained system 
is considerably more compressible than the corresponding atomistic one. (d) Despite recent efforts, the proper use of coarse- 
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grained potentials for dynamical simulations remains unclear. In particular the emergence of dissipation due to the co arse- 
grai ning step is mostly ignored, or, at best, included phenomenologically via some stochastic thermostat as done e.g. in lll23l 
Il35ll . For some notable exceptions see f58ll . Considerable work is definitely needed in order to arrive at a thermodynamically 
consistent description of a model system at the two levels of analysis (atomistic and coarse-grained), which will eliminate all 
these undesired symptoms and errors. 



B. Uncrossable chains of blobs 



Briels and collaborators 

CMEHIIMEU proposed a coarse-graining scheme wherein chains are sub-divided into a number 
of subchains of equal length; the center-of-mass of each such subchain is taken as the position of a corresponding mesoscopic 
particle called the blob. The blobs are connected by springs so that chain connectivity is preserved. Similarly to the coarse- 
graining procedure at the level of superatoms, the method makes use of a potential of mean force U mf = V(R^) for the 
position vectors of the n blobs, which ensures that the blob distributions in the atomistic and coarse-grained systems are the 
same. 

In order to describe shear flow effects in a velocity field of the from v x (r) = jr y , Kindt and Briels Il36ll proposed making 
use of the SLLOD algorithm [ 140]. Starting with a Langevin equation, such a method results in the following expression for the 
blob dynamics: 



Mi 



d 2 R t 
dt 2 



£Cff _ 



c 



ell 



dt 



C + 



(Ff • P, - jP lX P iy ) 



(2) 



where denotes the unit vector along the flow {x) direction, Mi is the mass of the z-th blob particle, R s ; its position vector, 
Ff = — dV/dHi the systematic force on particle i, £ the friction coefficient, and Ff the random force on particle i. Since the 
coarse-grained bonded and nonbonded interactions are so soft that unphysical crossing of two bonds woul d not be prohibited, 
equations (O are supplemented with an uncrossability constraint of the blob chains. Padding and Briels lll37l 1 1 4 ill realized 
this constraint by a method which explicitly detects entanglements and prevents chain crossings through a geometric procedure. 
The proced ure d etects possible chain crossings and defines an entanglement point X at the prospective crossing site. Padding 
and Briels Il36l Il37ll also proposed some non-trivial order- altering moves that lead to creation-removal of entanglements; 
these are important for the best possible realistic treatment of uncrossability constraints in simulations with the blob model. 
The Langevin equation of motion (O contains the blob friction coefficient £, whose calculation is not a straightforward issue 
even under equilibrium conditions. Despite this and its simplicity, the blob method has been found to capture correctly the 
viscoelastic properties of polymer melts with molecular length several times their entanglement length. From a numerical point 
of view, the method suffers from large requirements in CPU time, associated with the minimization algorithm for the location of 
entanglements which eventually limits simulations to chains made up of a finite number of blobs. 



C. Primitive paths 

A method to project atomistically detailed chains to smoother paths was proposed by Kroger et al. [84] through a projection 
operation that maps a set {r^}, i = 1, 2, . . . , N, of N atomistic coordinates of a linear discrete chain to a new set {Ri} of N 
coarse-grained ones defining a smoother path for the chain that avoids the kinks of the original chain but preserves somewhat 
its to pology (the main chain contour). The projection involves only a single parameter, £, whose value was obtained by Ottinger 
[142] by mapping the Porod-Kratky model (an atomistic model for a polymer chain) to a smoother chain with a Kuhn length 
equal to the entanglement length. In the limit of infinitely long chains, such a mapping suggests that the Kuhn length of the 
coarse-grained chain (the length of a segment between two entanglements) is equal to twice the tube diameter. The £-based 
method maps a particular chain onto a smoother path; however, the reduction of an ensemble of atomistic polymer chains to 
a mesh of primitive paths (PPs) as defined by the Doi-Edwards theory lfl4Tl requires that the projection satisfy not only chain 
continuity but also chain uncrossability. This subtle problem has been addressed only very recently through the seminal works 
of Everaers et al. ll82tl . Kroger ll8"oll . and Tzoumanekas-Theodorou |l43ll . Nevertheless, the simple £-based mapping has been 
very helpful in many respects; for example, it has allowed [85s, Il44ll to successfully calculate the zero shear rate viscosity of 
model polymer melts in the crossover regime from Rouse to entangled. 

The topological analysis of Everaers et al. 18211 is based on the idea that PPs can be identified simultaneously for all polymer 
molecules in a bulk system by: keeping chain ends fixed in space, disabling intra-chain excluded-volume interactions and retain- 
ing the inter-chain ones, and minimizing the energy of the system by slowly cooling down toward the zero Kelvin temperature. 
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This causes bond springs to reduce their length to zero, pulls chains taut, and results in a mesh of PPs consisting of straight 
segments of st rongly fluctuating length and more or less sharp turns at the entanglement points. The method can be modified 
ifsiL Il45l Il46ll to preserve self-entanglements or to distinguish between local self-knots and entanglements between different 
sections of the same chain. 

Kroger dSOfl also presented an algorithm which returns a shortest path and the related number of entanglements for a given 
configuration of a polymeric system in 2-D and 3-D space, based on geometric operations designed to minimize the contour 
length of the multiple disconnected path (i.e., the contour length summed over all individual PPs) simultaneously for all chains 
in the simulation cell. The number of entanglements is simply obtained from the shortest path as either the number of interior 
kinks, or from the average length of a line segment. Application of the algorithm to united-atom models of linear polyethylene 
(PE) 1 147] allowed the calculation of a number of important statistical properties characterizing its PP network at equilibrium 
and h elped make the connectio n with an analytic expression for the PP length of entangled polymers by Khaliullin and Schieber 
[148] following earlier works lll49[ Il50ll . A representative snapshot of the entanglement network as computed for a linear 
trans- 1,4-PB polymer (40 chains of C500 at T=450K and P=latm) with Kroger's method is shown in Fig. Ill Gl 

A third methodology for reducing chains to shortest paths has been presented by Tzoumanekas and Theodorou lll43ll where 
topological (chain uncrossability) constraints are defined as the nodes of an entanglement network. Through their contour re- 
duction topological analysis (CReTA) algorithm, an atomistic configuration of a model polymer sample is reduced to a network 
of corresponding PPs defined by a set of rectilinear segments (entanglement strands) coming together at nodal points (entangle- 
ments) by implementing random aligning string moves to polymer chains and hard core interactions. In addition to obtaining 
topological measures for a number of entangled polymers, Tzoumanekas and Theodorou II 14311 found that data for the normal- 
ized distribution of the reduced number of monomers (united-atoms or beads) in an entanglement strand are well described by 
a master curve suggesting a universal character for linear polymers. As analyzed by Tzoumanekas-Theodorou lll43ll . the master 
curve is also obtainable in terms of a renewal process generating entanglement events stochastically along the chain. 

Apart from some algorithmic details, the three methods lead to practically similar conclusions as far as the topological state 
of m any entangled linear polymer melts (PE, PB, PET and PS) is concerned. A significant result of all of them (see, e.g., 
lll43h is that the ensemble average of the number of monomers per entanglement strand Nes as computed directly from the 
topological analysis is significantly smaller than the corresponding quantity N c measured indirectly through N e = N (R 2 ) / (L) 2 
by assuming that PP conformations are random walks. This is due to directional correlations between entanglement strands 
along the same PP, which decay exponentially with entanglement strand separation. Therefore PPs are not random walks at the 
length scale of the network mesh size. 

D. Other single-chain simulation approaches to polymer melts: slip-link and dual slip-link models 

Doi and Edwards [ 151] and Doi and Takimoto lll52ll have proposed a description of an entangled polymer in terms of a slip- 
link model that can cumulatively account for chain confinement and constraint release in a consistent way. Slip-links do not 
represent an entanglement junction in real space; they are rather virtual links representing effective constraints whose statistical 
character is determined by other polymers. In the dual sli p -lin k version o f the model, the slip-link confines a pair of chains 
(and not a single chain). Masubuchi and collaborators l49l l50l Il53l Il54ll generalized the idea by regarding a slip-link as an 
actual link in real space. In their formulation, each polymer chain is represented as a linear sequence of entanglement strands 
considered as segments (phantom entropic springs) joining consecutive entanglement points (the beads) along the chain. These 
Rouse-like chains are all inter-connected by slip-links at the entanglement points to form a 3-D primitive chain network. The 
system is described by: the number Z of segments in each chain, the number n of monomers in each chain segment, and the 
position vectors of the slip-links or entanglement points. These state variables are postulated to obey certain Langevin- 
type governing equations in which the single relevant parameter of the primitive chain network is the average value (Z) of 
entanglements per chain. Defining the model functions and parameters on the basis of the results obtained from one of the 
three topological analysis discussed above leads to rheological predictions that follow quite satisfactorily experimental data for 
many polymer melts in shear but deviations are observed when the model is used to describe the elongational rheology of these 
systems. A generalization of the slip-link idea by Schieber and collaborators 11441 14611 to a full-chain slip-link model with a 
mean-field implementation of constraint release and constant chain friction (as opposed to constant entanglement friction) has 
been shown to provide accurate predictions of the G" and G" spectra for many polymer melts (such as PS, PB and PIB). 

E. Entire molecules 

In dilute polymer solutio ns, coarse graining polymer coil or star polymers to a a system of interacting soft particles has been 
explored in lfl55l fl56l SH . Kindt and Briels lll4ll Il58ll proposed such a type of coarse graining also for polymer melts where 
an entire chain is represented as a single particle. To account for the presence of entangleme nts w hich are considered to be 
responsible for the distinct viscoelastic properties exhibited by polymers, Kindt and Briels IU4UI158I1 introduced a second set of 
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variables, the number mj of entanglements between chains i and j. This governs the degree of interpenetration or overlapping of 
two chains whose centers-of-mass are fixed at a given distance. The state of the system is thus fully determined by the position 
vectors of the centers-of-mass of the A c h chains and the 7AT cn = A c h ( A c h — 1 ) / 2 entanglement numbers . The equilibrium 
density distribution function ^ for such a system is of the form: 



1 



exp 



fc R T 



^2-a(n tJ -n (ry)) 2 



(3) 



where k-& is Boltzmann's constant, T the temperature, U m f the potential of mean force, the double summation is over all 
interacting particle pairs, and a is a constant determining the strength of the fluctuations around a mean number no(ry) of 
entanglements between chains i and j. 7io(ry) is like an order parameter governing the "friction" felt by each chain and 
generating restoring elastic forces. According to Eq. (O, integration over Jiy results in a Boltzmann distribution for the A c h 
coordinates R;, thus the equilibrium statistics of the system is not al tered by the introduction of the entanglement s. Typica l 
expressions for U m f have been discussed by Padding and Briels lll36l Il37ll but also by Pagonabarraga and Frenkel Il59l Il60ll 
in their derivation of the "multi-particle dissipative particle dynamics" method. The method is capable of providing structural 
information only about the radial distribution function g(r) of the centers-of-mass of the chains. Representative results for a 
number of linear PE melts revealed a small correlation hole effect at the level of entire chains, which is consistent with data 
reported by Mavrantzas-Theodorou lll6lll through atomistic Monte Carlo simulations. No other signals of local structure could 
be discerned. Clearly, accounting for entanglements (which is necessary in order to produce the correct viscoelastic properties) 
seems to have a negligible effect on the structural properties at the level of entire chains. 

The single particle model has been proposed to describe systems where memory effects are dominant. This is the case 
for example of complex fluids involving colloid al particl es floating in a solvent in which a small amount of polymer is also 
dissolved. In this model, dynamics is described lll58lll62ll by generating (according to standard expressions for Smoluchowski- 
type equations) at every time step dt not only a displacement dRi in the position of each particle i but also a change dn^ in the 
number of entanglements. Memory effects are taken into account through transient forces: when two particles come together 
such that temporarily < no, their coronas are pushed apart causing a repulsion between the two particles. On the other hand, 
if the coronas are separated enough such that temporarily mj > no, the particles experience attractive forces. These phenomena 
cannot be studied by traditional Brownian Dynamics simulations where delta-correlated random displacements are assumed. 
Figure|4]shows the viscosity obtained from such a method referring to a typical resin with particles having a hard core diameter 
of 100 nm. The very same model has been successful in describing also shear banding and the chaining of dissolved colloids in 
viscoelastic systems lll63ll . 



F. Conformation tensor 



Based on the idea that, in a flow situation, certain structural variables depart from their values in the quiescent fluid while 
all other (faster) degrees of freedom are at equilibrium subject to the constraints imposed by the values of the slow variables at 
all tim es, coarse-gra ined models have also been developed where chains are described at the level of the conformation tensor c 
[|2l Il2l Il5l 1 1 64t Il65n . The latter might be defined via the average gyration tensor or via the tensorial product of the end-to-end 
vectors. The choice of c among the structural parameters marks a description in terms of a tensorial variable. In addition to a 
single conformation tensor one can envision a description in terms of many (higher-mode) conformation tensors, corresponding 
to the Rouse or bead-spring chain model. In general 016411 . from an A-mer chain, A(A — 1) /2 different conformation tensors 
Cij, i, j = 1, 2, . . . , N — 1 can be constructed, each one being identified as a properly dimensionalized average dyadic (QiQj) 
with Qi denoting the connector vector between mers i + 1 and i along the chain Q- 

For systems of entangled polymer chains, reptation theory suggests the distribution function of primitive path orientations 
as structural variable, see previous section. If chains contain long chain branches that can significantly affect the rheological 
response of the system, their contribution to overall system dynamics should also be accounted. This is the case, for example 
with H-shaped polymers for which a model built on two structural parameters, the tube orientation tensor S from one branch to 
the other branch along the chain, and a sca lar quantity A describing the length of the tube divided by the backbone tube length 
at equilibrium, has been proposed [24, 166]. This marks a description at the level of a tensorial and a scalar. 

The description in terms of a few, well-defined structural parameters (such as the conformation tensor, the configurational 
probability function, and the orientation tensor and a scalar) is very appealing because of the existence of well-founded models 
developed under the GENERIC framework of non-equilibrium thermodynamics II 1 6711 . In Section |IV] we will see that, among 
other things, this allows building thermodynamic ally guided multi-scale approaches by expanding the equilibrium statistical 
ensemble to incorporate terms involving conjugate variable(s) driving the corresponding structural parameters away from equi- 
librium. The GENERIC formalism here is an extremely useful tool, since it can guide us in linking the conjugate variable(s) 
to the applied flow field. It is not surprising therefore that this class of coarse-graining procedures constitutes one of the most 
understood and best-founded today. 
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G. Mesoscopic fluid volumes 



For the numerical simulation of flowing polymers, several mesoscopic models have been proposed in the last years that 
describe polymer (hydro-)dynamics on a mesoscopic scale of several micrometers, typically. Among these methods, we like 
to mention Dissipative Particle Dynamics (DPD) lll68ll . Stochastic Rotation Dynamics (sometimes also called multiparticle 
collision dynamics) |33il , and Lattice Boltzmann algorithms l30ll . Hybrid simulation schemes for polymer solutions have been 
deve loped recentl y, combining these methods for solvent dynamics with standard particle simulations of polymer beads (see 
70]). Extending the mesoscopic fluid models to non-ideal fluids including polymer melts is currently in progress 
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FIG. 2: (a) Polystyrene m and r diads in transplanar conformation (for clarity, hydrogen atoms on phenyl rings have been omitted), (b) 
Illustration of a mapping scheme fr om t he atomistic to a coarse-grained structure for PS wherein one bead corresponds to an m or r diad. 
[Reproduced with permission from 1 113], Figure 1.] 
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FIG. 3: (a) A snapshot of a fully equilibrated atomistic configuration of a 40-chain C500 trans- 1,4-PB melt at 500K and latm. (b) The 
corresponding entanglement network mesh as computed with Kroger's method [8Cj]. 



III. SYSTEMATIC AND THERMODYNAMICALLY CONSISTENT APPROACH TO COARSE GRAINING: GENERAL 

FORMULATION 



A. The need for and benefits of consistent coarse-graining schemes 



Under equilibrium conditions, statistical thermodynamics forms a bridge between thermodynamics (whose goal is the under- 
standing and prediction of macroscopic phenomena) and molecular physics (which focuses on the intermolecular interactions 
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FI G. 4: Viscosity versus shear rate for a typical resin, as obtained from the particle model of Kindt and Briels 1 158] and van der Noort ef 
al. Il62ll . The solid line represents experimental data while the crosses are simulation results including all forces in the stress tensor. The 
squares represent viscosities based on the conservative forces only. The circle at the vertical axis gives the zero-shear viscosity according 
to the Green-Kubo formula. In the inset, we show the overshoot in the instantaneous value of the viscosity at 7 =ls _1 . [Reproduced with 
permission from 1 162], Figure 3.] 



between the atoms making up the system). It provides therefore an interpretation of thermodynamic quantities from a molecular 
point of view. For a number of complex fluids (like colloids, liquid crystals, and polymers), information at a mesoscopic level 
of description (intermediate to molecular and macroscopic ones) is often extremely useful in understanding and predicting ma- 
terial behavior (see e.g. Sect. |II]and 1 12]). At equilibrium, statistical thermodynamics provides the framework for understanding 
system properties also at these intermediate scales. For example, as mentioned in Sect. Ill Al effective potentials describing 
interactions between coarse-grained (pseudo-)particles can systematically be derived by integrating out irrelevant degrees of 
freedom. 

Although extensions to capture dynamics at a mesoscopic level are in progress (see Sect. [TTJ>, most descriptions in terms of 
coarse-grained particles are so far largely restricted to equilibrium situations. For example, the recently proposed time rescaling 
approach (see IU14l,[l2"6[l ) for coarse-grained models does not seem appropriate to fully account for the increase in dissipation 
inherent in any meaningful coarse graining technique. We mention the case of hydrodynamic interactions in polymer solution 
which necessitate a description not in terms of a scalar frictional var iable but in terms of a tensorial friction matrix. 

Dissipation and friction are more properly accounted for in JH EM [HI [HI Wf% . These authors, however, arrive at a 
daunting assessment: "We therefore conclude that coarse-grained models lack thermodynamic consistency" Il38ll . As we will 
demonstrate below, for appropriately defined coarse-grained models, there is a way to restore thermodynamic consistency. So, 
contrary to the authors of 1 1 3811 . we believe that the recently introduced GENERIC formalism of nonequilibrium thermody- 
namics offers a framework for the development of true and complete coarse-graining strategies [11], in the sense that: (a) the 
resulting model is well-behaved and thermodynamic ally consistent, (b) it can be parameterized based on the information pro- 
vided by a lower resolution model, and (c) it can be improved based on microscopic simulations targeted to address the relevant 
structural variables and their dynamic evolution. Of course, a word of caution is in place here: coarse-grained models based on 
the GENERIC framework will rely on a number of strong assumptions (inherent to most projection operator based methods) 
implying a description in terms of a set of carefully chosen slowly-evolving state variables. The underlying assumption behind 
such a description is that of the existence of a clear time-scale separation between the evolution of these (slow) variables and that 
of the (eliminated) fast or irrelevant ones. The coarse-grained model of 013811 . for example, involves lumping ten beads along a 
chain into one or two blobs, for which the time-scale separation argument is questionable. Their negative conclusion about the 
thermodynamic consistency of the model is therefore not surprising. 
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B. Different levels of description and the choice of relevant variables 

Coarse-graining connects (at least) two descriptions of the same system at two different levels of detail: a low-resolution level 
and a higher-resoluti on le vel. We focus attention here to the case where the high-resolution level is the atomistic one, although 
this is not necessary Jl73ll . 

We consider a point in phase space z 6 V, where z = (ri, . . . , rjv; pi, . . . , Pat) is a short hand notation for the positions 
and momenta p^ of all N particles, at the microscopic level; this, for example, could be an all-atom or a united atom model or 
even the simpler and computationally more convenient FENE bead-spring model 11211 . All these three models are classified here 
as microscopic models due to the absence of dissipation and irreversibility. Dynamics at the microscopic model is governed by 
Hamilton's equation of motion 

i-J-£ (4) 
oz 

where H(z) is the microscopic Hamiltonian and J the symplectic matrix. Equivalently, Hamilton dynamics can be formulated 
by A = {A, H}, where {A, B} denotes the microscopic Poisson bracket between arbitrary functions A(z) and B(z). We 
recall that the basic properties of Poisson brackets are their anti-symmetry {^4, B} = —{B, A}, the L eibniz rule {AB, C} = 
A{B, C} + {A, C}B, and the Jacobi-identity, {A, {B, C}} + {B, {C, A}} + {C, {A, B}} = lfl67ll . 

Due to their very long relaxation times, there is a clear gap between the time scales that can be addressed in microscopic 
simulations of polymer melts and the relevant time scales in experimental studies. Although this prevents the direct applicability 
of brute-force microscopic simulations, it renders them ideal systems for a comprehensive understanding over multiple time and 
length scales by embodying the concept of multiscale modeling. The first and more important step in this context is the proper 
choice of the relevant variables at the coarser level. For simple fluids, densities of conserved quantities (mass, momentum, 
and energy) are the proper variables to consider if one is interested in hydrodynamic properties. For systems with broken 
symmetries, the corresponding order parameters constitute additional candidates for slow variables Il74ll . In the case of complex 
fluids, however, no general rules are available how the appropriate relevant variables should be chosen, and this emphasizes the 
importance of physical intuition 11561,1 10311 for the choice of variables beyond equilibrium. 

For polymers, one can be guided by available theoretical models. For example, orientational ordering in liquid crystals and 
liquid crystalline polymers can be described by the alignment tensor within the Landau-de Gennes theory Ill4ll . Birefringence and 
viscous properties in the case of unentangled polymer melts can be addressed by models based on the concept of a conformation 
tensor, see Sect. Ill Fl For branched polymers, a scalar variable is added to the conformation tensor in order to capture additional 
contributions to the stress tensor B2411 due to long arm relaxations. For entangled polymer melts, reptation theory provides a 
description in terms of a probability distribution function for the orientation of segments along the primitive path I12L Il4f1 . 

These theories are examples of mesoscopic or macroscopic models that lead to closed-form constitutive equations. Fur- 
thermore, they can all be described in the context of the single-generator bracket [ 175] or the GENERIC 1 167] formalisms of 
nonequilibrium thermodynamics, 

* = L.^ + M.£. (5, 
OX ox 

In ([5]), E(x) and S(x) are the coarse-grained energy and entropy functions, respectively. The anti-symmetric operator L defines 
a generalized Poisson bracket {A, B} = M • L • 4g| which shares the same properties as the classical Poisson bracket described 
above. The last term in (|5]l is new compared to Hamiltonian dynamics (HJl and describes dissipative, irreversible phenomena. 
The friction matrix M is symmetric 120211 and positive, semi-definite. Together with the degeneracy requirements L • (5S/8x.) = 
M • (SE/Sx) = 0, these properties ensure that the total energy E is preserved and S is not decreasing in time 016711 . 

The mesoscopic and macroscopic models come with a number of parameters, e.g. mean-field potentials, friction coef- 
ficients, effective relaxation times, etc., whose connection to molecular terms is not straightforward. It is the purpose of 
thermodynamically-guided, systematic coarse graining methods to address this issue. 



C. GENERIC framework of coarse graining 

Coarse graining implies a description in terms of a few, carefully chosen variables after the elimination of all irrelevant degrees 
of freedom . Inevitably, this comes together with entropy generation (irreversibility) and dissipation. In the GENERIC formalism 
Jl67lll73ll . the emphasis is then shifted from the fundamental time evolution equation itself to the individual building blocks of 
that theory ([5]), and opens up the way toward the development of consistent coarse-graining strategies. T he interest ed reader is 
referred here to Ref. [11]. More details on the statistical mecha nics of coa rse graining can be found in lf56lll67lll76ll : for coarse 
graining of simple fluids within the GENERIC framework, see Il77lll78ll . 



11 



1. Mapping to relevant variables and reversible dynamics 

For every microstate z of the system, the instantaneous values of the relevant variables are defined by a set of phase space 
functions II (z). The functions II(z) cannot generally be identified with x; they are rather connected with x through x = {n(z)), 
i.e., as averages based on a suitable probability density p x (z) at the microscopic phase space T. Thus, the coarse-grained energy 
-E(x) is obtained from the microscopic Hamiltonian H(z) by straightforward averaging, 

B(x) = (H(z)) x , (6) 

and, similarly, the coarse-grained Poisson bracket is obtained from the average of the classical Poisson bracket, 

{A,B} = ^-L kl -^-; L H (x) = ({n fc ,n I }) !B . (7) 

OXk OXl 

Eqs. (O and (0 define the reversible part of GENERIC (O in terms of a coarse-grained Poisson bracket lll74ll79ll . The addi- 
tional terms related to dissipation and increase in entropy have to be accounted for by the irreversible contribution to GENERIC 
(0 and are described just below. 



2. Irreversibility and dissipation through coarse graining 

The fact that we do not account explicitly for the irrelevant variables at the level of the GENERIC framework leads to entropy 
increase and additional dissipation at the coarser level of description 111 ill . 




FIG. 5: (a) Schematic illustration of freely jointed chain with end-to-end vector R cc . (b) Fluctuations of polymer chain in shear flow around 
stationary state. Ellipses indicate the eigenvalues and orientation of eigenvector of (R cc R co ). [Figure courtesy of M. Kroger, ETH Zurich.] 

To illustrate this, let us consider the simple example of a freely jointed chain, shown schematically in Fig.|5^. We can describe 
such a chain by the set of all connector vectors Qj = Tj+i — rj, j = 1, . . . , N. All admissible configurations with \Qj\ = b 
have equal probability. If we decide choose the end-to-end vector R oc as the only relevant variable, then there are in general 
many configurations {Qj} which are compatible with R cc . The coarse-grained entropy S'(Ree) is a measure of the number of 
these configurations. One finds that the probability of R cc is Gaussian around R cc = 0, which implies that there are many more 
coiled configurations compared to stretched ones. The associated entropy S^Roo) = S(0) — 3H^ c /2Nb 2 decreases for chains 
undergoing stretching and leads to a restoring force which is known as "entropic spring" in coarse-grained polymer models 01 . 
This illustrates the emergence of (additional) entropy through coarse graining. 

We turn now to the discussion of the probability distribution p x (z). In sharp contrast to equilibrium statistical mechanics, 
there are unfortunately no general results for the probability distribution of nonequilibrium states. Even fo r non equilibrium 
stationary states there are at present only a few results for very special model systems available (see e.g. lll80l0 . Systems, 
however, where the time-scale separation assumption holds are well described within the quasi-equilib rium appro ximation that 
treats the nonequilibrium system as an equilibrium one for the present values of relevant variables ill 8ll Il82l Il83ll . In the 
generalized microcanonical ensemble, all microstates z that are compatible with given values of the relevant variables II (z) have 
equal probability. The corresponding entropy is a measure of the number of such microstates z that are compatible with a given 
coarse-grained state. For practical calculations, it is more convenient to pass to the generalized canonical distribution. In analogy 
to equilibrium statistical mechanics, the average values Xk = (Hk{z))x are constrained to prescribed values with the help of 
Lagrange multipliers The generalized canonical distribution can then be obtained from the maximum entropy principle and 
reads 

e -£ fc A fe n fc (z) 

Px(z) = J- rdge -S.*n,(.) < 8 > 
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where the A^'s have to be chosen so as to satisfy Xk — Jdz Uk(z)p x (z). The quasi-equilibrium entropy associated with Eq. ([8]) 
is 

5(x) = fc B V A fc x fc + k B In f dz e- £' A,n, («) (9) 

The coarse-grained entropy plays the role of an effective potential for the relevant variables. Determining the functional form of 
S(x) from © presents a challenge, since the explicit expression for the Lagrange multipliers A^(x) is in general unknown. A 
successful method for extracting at least partial information on S(x) has been explored in Il64ll from atomistic simulations of a 
polymer melt in elongational flow. We discuss this issue further in Sects. IrVBl and lTVCl 

Having specified the nonequilibrium ensemble and coarse-grained entropy, we finally like to discuss the increase of dissipation 
through coarse graining in more general terms. We have seen above that many microstates (values of connector vectors) are in 
general compatible with a given coarse-grained state (defined by the value of the end-to-end vector in the above example). 
Conversely, this implies that a coarse-grained state does not uniquely determine the microstate. The dynamics on the coarse- 
grained level has necessarily a stochastic character which is known as fluctuations, see Fig.|5p. If those fluctuations are correlated 
in time, they have to be accompanied by dissipation, as required by the fluctuation-dissipation theo rem |[56ll . These qualitative 
observations are put into a solid theoretical framework by the projection operator formalism lll03ll . It should be emphasized 
that projection operators provide exact relations for any set of variables in terms of complicated integro-differential equations. 
Simpler, closed-form equations without a memory integral, however, r esult only in cases when the time-scale of the chosen 
variables is well-separated from those of the irrelevant ones Il03[ fl841. A prominent example where this assumption seems 
not to be met is the dynamics of glassy polymers, where usually mode-coupling approximations for the memory kernel are 
employed 1 1 85l 1 1 86ll . We here insist on the time-scale separation, which severely restricts possible choices of relevant variables 
where such a separation can hold. For glassy polymers and glasses, in general, an approp riate set of relevant variables is not 
known at present, although some promising first steps have been taken recently lll87ll 18811 . These restrictions are the price to 
pay for a proper coarse-grained description with a well-defined entropy and without accounting for memory effects. In this case, 
the dissipation matrix M as derived from the projection operator formalism reads 

M kl = -L f S di(ri f fc (i)n[(0)), (10) 

where ri^ is the fast part of the time derivative of the macroscopic variables lll67[ Il73l Il82ll . The separating time scale r s 
should be chosen large enough to comprise all the fast fluctuations that are not captured on the coarse-grained level Il67ll . Thus, 
the friction matrix M arises due to fast fluctuations that are not resolved at the coarser level. The numerical evaluation of the 
dissipation matrix (TToT > for a polymer melt is described in Sect. |IV C| 



IV. THERMODYNAMICALLY GUIDED COARSE-GRAINED POLYMER SIMULATIONS BEYOND EQUILIBRIUM 

A. GENERIC coarse-graining applied to unentangled melts: foundations 

Unentangled polymer melts are usually described in terms of the conformation tensor which provides an overall picture of the 
entire polymer chain. From the point of view of the GENERIC formalism, this implies a description where, in addition to the 
hydrodynamic fields mass p, momentum g, and energy density e, the conformation tensor c is also included in the vector of state 
variables, 

x=(p,g,e,c). (11) 

Let z = (ri, . . . , tn; Pi, ■ ■ ■ , pjv) denote the microstate defined by the positions and momenta of all particles. The chosen 
macroscopic variables x = (n(z)) are defined as follows. The mass density is defined by 

p(r;t) = (^m^(r -!•,■(*))) = (H p ), (12) 

3 

where rrij denotes the mass of particle j. Similarly, the momentum density is obtained by 

g(r;i) = Epjjfr-rjft))) = (II 9 >. (13) 

j 

From p and g, the macroscopic velocity field v(r) is defined by v(r) = g(r) / p(r). The total energy density can be expressed as 

e(r;t) = e^(r - r,(i))> = <II e ), (14) 

3 
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where e. 



c 3 — (l/2)mjUj+&j with Uj = pj/nij — v(rj) the peculiar velocity of particle j and the potential energy of particle 
j. Finally, the additional, internal variable c is a symmetric, second-rank tensor which is defined by 



1 ^ch 

c ^ *) = — E^(r - r a c (t))) = (n c >, 

a— 1 



(15) 



where iV c h is the number of chains in the system and N a = me number of particles in chain a. The center of mass 

of polymer a is denoted by r° = N~ z^j G / rj. The tensor II is a conformation tensor of a single chain and quantifies the 

instantaneous, internal structure of polymer a. Examples are the gyration tensor II = iV^ 1 J2jei ( r j — r c)( r j — r c) or me 
tensor product formed either by the end-to-end vector or the first Rouse mode. 

The macroscopic energy E is obtained by straightforward averaging of the microscopic Hamiltonian, see Eq. (0. 

The resulting expression for the distribution function of the generalized canonical ensemble reads 



p x (z) = Z" 1 exp [-/3pV - - A: II C ] 



(16) 



where (3 = (/cbT) -1 , V the volume o ccupi ed by the N pa rticles, $ the total potential energy, Z the normalization integral and 
p the pressure (see e.g. Sect. 8.2.3 in lll67ll and lll64 Il89ll where ( TTol l is used in nonequilibrium situations). The macroscopic 
entropy associated with the generalized canonical distribution is given by Eq. ©, which here reads 



S(x) = S (T,V,N) + kt 



In 



V N J 



+ /3pV + /3{<P) x + \: 



(17) 



where V — (V) x is the average volume and So the entropy of an ideal gas of N particles. In addition to the usual Lagrange 
multipliers (3 and f3p that are associated with total energy and the volume (for homogeneous density), respectively, the additional 
Lagrange multiplier A is identified as 



OC 



(18) 



For a numerical calculation of the Lagrange multiplier A for a model polymer melt see Sect. IIV Bl and llV Cl 

The matrix L defining the coarse-grained Poisson bracket (0 is obtained by inserting the defini tions (112l|15t of the coarse- 
grained variables into Eq. (0. Details of the straightforward calculations are presented in Ref. IU90I1 . 

From the degeneracy requirement on the Poisson bracket {S, E} — mentioned in Sect. IIIIB1 one finds that the entropic part 
of the macroscopic stress tensor has to be of the form 



ds 

cr = -Pcffl - 2Tc • — , 

OC 



(19) 



where jp e ff is the effective scalar pressure and s the local entropy density. The same form ( fT9] l has been previously found in 
As far as the dissipative bracket and the associated friction matrix M are concerned, a direct calculation of the fast time 

• f 

evolution IF appearing in Eq. ( fTOb shows that 11^ = 0, II g = V • <r tot , where <x tot is the instantaneous value of the total 
stress tensor. The expression for 11^ containing the heat flux and viscous heating can be found in [ 190]. The integral of the time 
correlation function of these fast fluctuations that appears in Eq. (fTOb can in most cases only be determined numerically. How 
these quantities can be extracted from molecular dynamics simulations for a model poly mer melt is des cribed in Sect. IIV CI 
The resulting GENERIC equations (0 for the present choice of relevant variables are lll67lll90lll9lll 



dl P 



dt 9a 



*/3,[l] 



de 

dc 



-T 



ds 



dCau 



(20) 



where cm denotes the upper-convected derivative of c, cm = <9 t c + vc — ck t — kc, k = (Vv) T the transpose of the velocity 
gradient, and e and s the ener gy and ent ropy density, respectively. Since we consider in the following only isothermal conditions, 
the reader is referred to Refs. Tl90l OH for the rather lengthy expression of the internal energy balance. Furthermore, additional 
second order dissipative processes appearing in Eq. d20i > are discussed in Ref. 1 190]. 
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The macroscopic stress tensor appearing in the momentum balance equation (l20b is given by 

-Ic^k, (21) 

where C ( - <T<T ' = f^dt (<x f (t)cr f (0)) x is a Green-Kubo formula for the viscosity contribution of fast (on time scale shorter than 
t s ) stress fluctuations. This finding is in agreement with previous simulation studies on b ead-spring chain polymer melts that 
found it necessary to include a simple fluid background viscosity in their analysis lf57lll92ll . 



-p e gl + 2c ■ 



dc dc 



B. Thermodynamically guided atomistic Monte Carlo methodology for generating realistic shear flows 

We discuss here how one, guided by principles of nonequilibrium thermodynamics, can make use of the Monte Carlo technique 
to drive an ensemble of system configuration s to sample statistica l ly a ppropriate steady-state nonequilibrium phase-space points 
corresponding to an imposed external field Il6ll Il64l Il93l Il94l Il95ll . For simplicity, we limit our discussion to the case of 
an unentangled polymer melt. The starting point is the probability density function p x of the generalized canonical GENERIC 
ensemble (fTol i for the same set of slow variables (flTT i as in Section llV Al Then, following Mavrantzas and Theodorou lll6lll . we 
extend the Helmholtz free energy, A, of equilibrium systems to nonequilibrium systems as 

= ~dT + M d - k B TX: dc (22) 

where fi is the chemical potential. The last term accommodates the effect of the external field (e.g., a flow) for which A represents 
a nonequilibrium force variable conjugate to c. According to Eqs. (TToT l and (l22l . one can carry out Monte Carlo simulations in the 
expanded {N^nPT X} ensemble exactly as in the corresponding {N^nPT} equilibrium ensemble (with n denoting the total 
number of atoms in the system) by assigning non-zero values to the field A. This is the key point of the new method opening up 
the way toward sampling steady-state nonequilibrium phase points of the system corresponding to a given flow field with Monte 
Carlo by suitably choosing the components of A. For the case of a simple shear flow, for example, from the symmetry property 
of c, we recognize that A is to have only four independent non-zero components: X xx , X xy , X yy , and X zz . In order to specify 
their numerical values for a given shear rate 7 (these are needed to be used as input in the GENERIC MC simulations), one can 
resort to the fundamental GENERIC evolution law for the set of state variables x. Based on this and Eq. ( TToT l for the definition of 
A, we see that, indeed, we can assign a kinematic interpretation to the Lagrange multiplier A, since for a nonequilibrium system 
that has reached a steady state, Eq. @ simplifies to 

A=-Im-L(x).^ (23) 
«b ox 

For example, for all known conformation tensor viscoelastic models, the corresponding evolution equation for the conformation 
tensor reads: 

SA(c) 1 SA(c) 
c a a,[i] = -A a( 3 7 e-^ = -nfc B TA a/ 3 7£ Q; 7e ; <x a p = — — . (24) 

OC-ye nKBl 0C a f3 

where (for simplicity) we have replaced the tensor A with the tensor a defined through A = —(N/V)a. In Eq. ( f24b . c a p 
denotes the upper-convected derivative of c a p and n the chain number density, and the Einstein summation convention has been 
employed for repeated indices. Note also that, in the case considered here, the element M44 of the M matrix in Eq. d23l has the 
form of Tk a p ie (see Refs. II 1 6411 and II 19311 for details) where the fourth-order relaxation matrix A for most single-conformation 
tensor models can be cast into the following general form: 

A a/ 3 7£ (c) = h{h) ( as T Cp £ O a ^ 

) + 2/a(Ji) ( c/a 7 ) (25) 

where I\ is the first invariant of c (i.e., the trace of c), 5 the unit tensor, and fx and fi arbitrary functions of l\. With the help of 
Eqs. (l24l and (|25l l, for the case of a steady-state flow described by the kinematics 



Vv= I 7 I (26) 



we find that the form of a that generates shear is 







<-xy 

*yy 



(27) 
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Although, however, nonequilibrium thermodynamics has helped us define the functional form of a, the exact relationship 
between its three non-zero components (a xx , a xy , and a yy ) on the applied shear rate 7 remains still undetermined. One way 
to come around this problem is to make explicit use of specific expressions for the matrix M proposed by GENERIC for a 
viscoelastic model. In such a case, however, the results will be model dependent and not representative of the true structure 
developed in the system in response to the applied shear rate 7. Baig and Mavrantzas II 19311 proposed overcoming this by 
computing a iteratively so that, for a given value of 7, the resulting average conformation of the simulated melt is the same as 
that predicted through a brute force application of the NEMD method. 

Baig and Mavrantzas 119311 demonstrated the applicability of such a hybrid GENERIC MC - NEMD approach for a relatively 
short unentangled PE system, C50H102, for different nonequilibrium states corresponding to different values of the Deborah 
number (De). De is defined as the product of the imposed shear rate 7 and the longest relaxation (Rouse) time of the system, tr, 
at the temperature and pressure of the simulation. If x is the flow direction and y and z the velocity gradient and neutral directions, 
respectively, then the three non-zero components of a can be computed iteratively so that the values of the conformation tensor 
c for the system (at the given value of De) from the GENERIC MC and the NEMD methods coincide. Representative results are 
shown in Figures[6]and[7] Figure |6]presents the values of the non-zero components of a. that were found to reproduce accurately 
the corresponding nonequilibrium state for the simulated C50H102 system as a function of the imposed De. Figure|7J on the other 
hand, presents comparisons of the conformation tensor between the GENERIC MC simulations (corresponding to the a-values 
shown in Figure|6]l and the direct NEMD simulations, confirming that c xx , c xy , and c yy from the GENERIC MC and the NEMD 
simulations, respectively, superimpose. It is only for the c zz component of the conformation tensor that Figure [7] reveals an 
inconsistency between the two methods. As argued by Baig-Mavrantzas [ 193], this is related to the selection of a zero value for 
the a zz component of a., as suggested by the general expression, Eq. (f25t . To reproduce exactly also the zz-component of c, 
a non-zero a zz component should be incorporated in the GENERIC MC simulations. This is a significant accomplishment of 
the new methodology, since it suggests that the rather general form of the friction matrix, Eq. d25l l, for this conformation tensor 
family of models is not complete. As demonstrated by Baig-Mavrantzas in a recent publication 1 196], this can be achieved by 
including in the relaxation matrix a terms beyond the symmetries implied by Eq. ( 1251 ). without violating the Onsager-Casimir 
reciprocity relationships nor the 2nd law of thermodynamics. 
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FIG. 6: Plot of the thermodynamic force field, a, vs. De number for the C50H102 PE melt (T=450K, P=latm). [Reproduced with permission 
from [1930 , Figure 1]. 

The information provided by the GENERIC MC simulations is important in many aspects: 

• The dependence of the components of the tensor a. on De is directly related to the (nonequilibrium) free energy of system 
- see Eq. (1221 1. Therefore, with the proposed methodology one can accurately calculate the free energy of the simulated 
system by requiring a series of simulations, according to the thermodynamic state points, by varying one component of 
a and fixing the rest and then using thermodynamic integration. This can serve as a starting point for developing more 
accurate viscoelastic models. 

• The new thermodynamically-guided method can help overcome the problem of long relaxation times (and of statistical 
noise in the Newtonian plateau) faced in brute-force NEMD simulations by providing the initial configuration at the 
relevant nonequilibrium state for a given De. 

• The new method can also be combined with recently proposed coarse-graining simulation strategies for long polymer 
melts to enable the simulation of the viscoelastic properties of high molecular weight polymers, comparable to those 
encountered in practical polymer processing. 
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FIG. 7: Comparison of the conformation tensor c components between NEMD and GENERIC MC simulations, as a function of De number: 
(a) c xx , (b) c xy , (c) c yy , and (d) c zz . The error bars are smaller than the size of the symbols. [Reproduced with permission from ll93ll . Figure 
1]. 



C. Systematic time-scale bridging molecular dynamics for flowing polymer melts 

We consider again a description of the polymer melt coarse-grained to the level of the conformation tensor. The corresponding 
Poisson bracket is known analytically, see Sect. lIV Al Same as in Sect. lIV Bl we investigate the nonequilibrium stationary state 
of the polymer melt in a given flow situation, and therefore face the same problem of solving the stationary GENERIC equation 
(|23| > self-consistently. Here, we complete the studies reported in Sect. lIV Bl and consistently determine also the friction matrix 
M from microscopic fluctuations according to general formula ( fTOb . Our presentation mainly follows Ref. 1 197]. 
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FIG. 8: Schematic illustration of systematic time-scale bridging algorithm which consistently combines Monte Carlo and Molecular Dynamics 
simulations. 

1. Systematic time-scale bridging algorithm 

The coarse-grained energy and entropy, as well as the Poisson bracket require only static information and can therefore be de- 
termined very efficiently by Monte Carlo (MC) simulation methods, see Section llV Bl and Il8lll98ll . Only the friction matrix M 
depends on dynamical properties, thus its numerical evaluation requires dynamical simulation, in our case Molecular Dynamics 
(MD). The GENERIC coarse-graining approach therefore suggests to combine the strengths of MC and MD simulations in a 
well-defined way to break the time-scale gap between microscopic and macroscopic scales. In order to implement these ideas 
consistently we propose a hybrid algorithm [ 197] schematically illustrated in Fig.[8]as a general strategy for time-scale bridging 
simulations based on G ENERIC. For the special case of flowing, unentangled polymer melts, the algorithm was implemented 
and tested in Ref. lll97ll . For simplicity and speed of calculations, the classical FENE bead-spring model introduced in l60ll was 
used in these studies, although making use of an atomistic model does not pose any extra difficulties. 

For this model system subject to a stationary flow with fixed velocity gradient k, the algorithm illustrated in Fig. [8] can be 
implemented as follows. 

1. Start with an equilibrium system with k = A = 0. 

2. Use a Monte Carlo scheme in order to generate an ensemble of n s (typically n s — 500) independent configurations that 
are distributed according to the generalized canonical distribution (IT6b with the current value of A. Calculate the value of 
the relevant variables in this ensemble, x = (II)^. In order to effi cient ly generate an ensem ble of n s such configurations, 
a slight modification of the Monte Carlo algorithm proposed in IU99I1 was used in II 9711 . The numerical values of the 
coarse-grained variables x = (n}^ can then be estimated as the ensemble average x = (l/n s ) Y^kLi n(z^) of the n s 
configurations {z/j}. In the last step, Maxwellian distributed velocities are assigned to the particles, realizing equilibrium 
in momentum space for the present choice of (velocity-independent) relevant variables. 

3. The Monte Carlo generated ensemble is used as initial condition for Molecular Dynamics simulations of Hamilton's 
microscopic dynamics. We use a standard velocity-Verlet algorithm that preserves the symplectic structure to simulate 
trajectories Zfe(t), k = l,...,n s , during a "short" time interval < t < t s . The separating time scale t s is short 
enough, such that the relevant variables x do not change significantly during the MD simulation. For this reason, the MD 
part of the simulation does not need any constraints such as thermo- or barostats nor flow-adapted boundary conditions. 
Performing short time, unconstrained, microcanonical molecular dynamics simulations is one of the great benefits of the 
present approach as it makes the scheme both highly efficient and applicable to arbitrary flow situations that - due to the 
lack of corresponding boundary conditions - could not be simulated so far. 

4. From the particle trajectories Zfc(t), we evaluate the friction matrix M from Eq. ( fTOb . We make use of time-translational 
invariance to equivalently rewrite Eq. ( fTOb as 

M = (M(z)) x , M(z) - -^A Ts n(z)A Ts n(z) (28) 

2k B T s 

where A Ts .n(z) = II(z(t s )) — n(z(0)) denotes fast fluctuations of II (on the time scale r s ). Equation d28l is more 
convenient for numerical evaluation than HOi . 

5. Updated value of the Lagrange multiplier A are calculated from the stationary GENERIC equation (1231 by inverting the 
symmetric, positive semi-definite matrix M. 
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FIG. 9: The friction matrix calculated from (A) : Eg. 4 10b compared to (B): the values obtained from Eq. l |28t , as a function of the separating 
time scale r s . [Reprinted with permission from 1 197], Figure 2.] 

6. The procedure is now repeated until consistent values x, M, A for given k are obtained. Alternatively, one may use an 
efficient reweighting scheme if k is changed only slightly and A is already close to the true value A — > A + S\. Then, the 
explicit form of the generalized canonical distribution can be exploited to solve the nonlinear system of equations 

" , A ll ;.. 

= 2^ [Rfc + k B M{z k ):5\] w k , w k = — _ sx . n(z ; (29) 

fc=l ^*' e ' k " 

for S\. This first order scheme solves the stationary GENERIC equations ((23), where = k ■ n(zfc) + n(z&) • k t + 
fcBA^(zfc) : A is the error in the previous value of A. In a shear flow, for example, Eq. ( |29l represents six equations 
and six unknowns. The solution <5A of (|29l allows one to calculate the reweighted slow variables and friction matrix, 
x = J^k u 'fcn(zfc), M = J2k u 'fcAl(zfc), as well as updated Lagrange multipliers, A — > A + <5A. Finally, the flow rate k 
is increased, and the procedure is started again, until the control parameter space has been swept through. 

With such a scheme, we establish the coarse-grained model along one-dimensional paths in the parameter space. Choosing, 
for example, viscometric flows of varying strength k is analogous to the situation encountered in experiments. 

2. Fluctuations, separating time scale, and friction matrix 

We have already emphasized several times that "fast" but correlated fluctuations give rise to dissipation on the coarse-grained 
level of description, which is described here by the friction matrix M, Eq. (TTOb or (l28l i. The notion "fast" is defined here 
by times t smaller than the time scale r s , that separates the evolution of the relevant variables x from rapid dynamics of the 
remaining degrees of freedom. The existence of such a time scale (which is equivalent to the crucial assumption of time-scale 
separation discussed in Sect. UTIb is not obvious. Here, we observe that the correlation functions Cj-iif) — (IT|(t)n[(0)) x decay 
monotonically over a few molecular (Lennard- Jones) time units r. This shows that those fast fluctuations are indeed correlated 
only over short times compared to typical polymer relaxation times (which are huge relative to r). Therefore, we find that 
the friction matrix, which is proportional to the integral over C(t), rapidly converges towards a value that is approximately 
independent of r s in a broad range 5 < r/r s < 50, see Fig. [9] 

3. Results 

Before discussing the results obtained with the proposed time-scale bridging algorithm, we like to mention several consistency 
checks that can be performed in order to test the range of applicability of the coarse-grained model. First, we compare two 
expressions for the macroscopic stress tensor. One is the standard virial expression er = — V^ 1 (rF) x , where r and F are 
the relative position and forces between particles. The kinetic contribution is found to be negligible in dense systems such as 
polymer melts as long as the flow rates are not too high 111 211 . Evaluating the expression for the stress tensor er in the generalized 
canonical ensemble leads to the expression er p = — 2l / ~ 1 fceTx • A for the (entropic part of the) polymer contribution to er, see 
Eq. ( fT9l . From Sect. |IV A| we know that er and er p differ by a simple-fluid contribution. Accounting for this contribution via 
the non-bonded short range repulsive interactions, we have verified that the two expressions for the stress tensor agree with each 
other for the flow rates studied. Next, by its definition and the symmetry of II, the matrix M possesses some basic symmetries 
that can be used to test the statistical accuracy of the ensemble averages. Finally, for the case of simple shear, k — -ye^e,,, we 



19 




-4 -35 -3 -25 -2 -IS 

1="3 1C| 5 hear rate 



FIG. 10: Polymer contribution to the shear viscosity as a function of shear rate for different molecular weights. Reference results 1 60] obtained 
with standard NEMD simulations for N — 30 are indicated. Inset: Zero shear rate viscosity r/o and first normal stress coefficient as a 
function of chain length. The expected scaling r/o oc N and ^i oc N 3 is observed, as shown by the dashed lines. [Reprinted with permission 
from 1 197], Figure 5.] 

have used the identity (xn — X2i)x^ = (An — A22) A^ 1 , that can be derived from the stationary GENERIC equations lll93ll . in 
order to check the consistency of our results. In our studies, this identity holds within error margins for the flow rates considered. 
We observed that breakdown of this relation at high flow rates signaled problems with the coarse-grained model as it can no 
longer capture the relevant dynamical processes at these elevated rates. 

For the case of simple shear flow, we validated the algorithm by reproducing the chain-length dependence of the zero shear 
rate viscosity and of the first normal stress coefficient, which are known in the literature i60ll . see Fig.[l0] Also the shear rate 
dependence of the v iscosity obtained with the time-scale bridging algorithm is in very good agreement with standard NEMD 
results [60 s l200l l20lll . as shown in Fig. [10] More results can be found in Ref. lll97ll . As mentioned above, the flexibility of 
our time-scale bridging simulations allows us to study arbitrary flow fields. We therefore could perform the first steady state 
equibiaxial simulation for polymer melts. Results for this as well as other elongational flows can also be found in lll97ll . 

V. CONCLUSIONS AND PERSPECTIVES 

The tremendous multiplicity of length and time scales in polymeric systems clearly calls for systematic, multiscale modeling 
approaches in which a higher-resolution model is consistently coupled with a lower-resolution one. In particular, if one is 
interested in describing relaxation processes and structure development under nonequilibrium conditions, most present-day 
coarse-graining strategies based on the use of effective potentials are of limited use since they do not account for the additional 
dissipation and irreversibility accompanying inevitably the elimination of fast degrees of freedom in favor of a smaller set of 
slowly-relaxing structural variables. Therefore, thermodynamically guided simulations are very important and useful, where 
one takes full advantage of the underlying principles of nonequilibrium thermodynamics and statistical mechanics. There, the 
emphasis is shifted from the time evolution equations (which respect important physical laws such as the Onsager reciprocity 
relationships for the transport coefficients and the 2nd law of thermodynamics) to its four building blocks: the energy E, 
the entropy S, the Poisson matrix L and the friction matrix M, describing the reversible and dissipative contributions to the 
dynamics. 

We have outlined such a methodology for the case of unentangled polymer melts for which, guided by network theory ap- 
proaches to polymer elasticity, the appropriate coarse-grained variable x is the conformation tensor. The underlying, microscopic 
model is simulated by the nonequilibrium molecular dynamics (NEMD) method. The relevant nonequilibrium state is assumed 
to be given by a generalized canonical distribution incorporating a conjugate variable (the Lagrange multiplier) A to the con- 
formation tensor. Monte Carlo simulations in this ensemble are employed in order to calculate the values of the slow variables 
x and the static building blocks E and L. For a given value of imposed flow rates, the Lagrange multiplier can be determined 
iteratively so that the solutions of the micro- and macro-solvers for the coarse-grained structural variables coincide. Through 
this one can compute model-independent values of the Lagrange multiplier, which for a wide range of strain rates (covering 
both the linear and the nonlinear viscoelastic regimes) bring results for the overall polymer conformation from the two models 
(microscale and macroscale) on top of each other. We presented two approaches to obtain the missing blocks of the macroscopic 
model. In Sect. lIV Bl the computed values of the Lagrange multiplier are compared with those corresponding to specific choices 
of the friction or relaxation matrix M in the macroscopic GENERIC model (addressing the chosen structural variable; here the 
conformation tensor), one can identify shortcomings and suggest improvements. And this is the biggest advantage of the new 
framework since the multiscale model proceeds without a priori knowledge of the exact form of the macroscopic model. Being 
built on the GENERIC framework of nonequilibrium thermodynamics, what is only needed is just to rely on the nature of the 



20 



chosen structural variables at the coarse level. This, further, emphasizes the significance of the choice of variables in the method. 
In the second approach, Sect. |IV CI we introduced a novel, low-noise, time-scale bridging strategy for the same system (low 
molecular weight, unentangled polymers) subjected to homogeneous flow fields. Through an alternating Monte Carlo-molecular 
dynamics iteration scheme we were able to obtain the model equations for the slow variables. For a chosen flow (including elon- 
gational ones), the method predicts structural as well as material functions beyond the regime of linear response. The method 
is simple to implement and allows for the calculation of time-dependent behavior through quantities readily available from the 
nonequilibrium steady states. In the end, it is only when all three different methodologies (macro-model, micro-model, and 
the macroscopic viscoelastic GENERIC equation bridging them) come together to complement each other that the entire multi- 
scale strategy can be considered as successful. Then, simulation techniques are elevated from brute-force computational tools 
to sophisticated techniques capable of mapping the detailed description of the system to a handful of carefully chosen variables 
whose dynamics (time evolution) is also faithfully described by an accurate analytical model. 

Future efforts will address other systems such as entangled (linear and branched) polymers where, inspired by the correspond- 
ing GENERIC formalism, one should resort to a description in terms of the orientational distribution function of an entanglement 
segment along the primitive path of the chain. 
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